Logistic and Autologistic in R

library(sf)
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(RColorBrewer)

sl <- read_sf("./SriLankaNYCFinal/SriLankaCorrectedFinal_UTM.shp")

#create a poor Secretariat on not variable
sl$POOR[is.na(sl$POOR)] <- 0
sl$POORFac <- factor(sl$POOR)
#ESDA
#map poor or not
plot(sl["POOR"],main='Poverty by Secretariat (binary)', breaks="jenks")

#map percentage poor
plot(sl["PCTPOOHH"],main='Poverty by Secretariat', breaks="jenks")

#histograms to see the distribution
hist(sl$POOR)

hist(sl$PCTPOOHH, breaks=20)

summary(sl)
##     OBJECTID          ID_0         ISO               NAME_0         
##  Min.   :36745   Min.   :217   Length:236         Length:236        
##  1st Qu.:36822   1st Qu.:217   Class :character   Class :character  
##  Median :36898   Median :217   Mode  :character   Mode  :character  
##  Mean   :36895   Mean   :217                                        
##  3rd Qu.:36967   3rd Qu.:217                                        
##  Max.   :37033   Max.   :217                                        
##       ID_1          NAME_1               ID_2           NAME_2         
##  Min.   : 2.00   Length:236         Min.   : 20.00   Length:236        
##  1st Qu.: 6.00   Class :character   1st Qu.: 96.75   Class :character  
##  Median :12.00   Mode  :character   Median :172.50   Mode  :character  
##  Mean   :12.13                      Mean   :169.53                     
##  3rd Qu.:17.00                      3rd Qu.:242.25                     
##  Max.   :23.00                      Max.   :308.00                     
##     HASC_2              CCN_2      CCA_2              TYPE_2         
##  Length:236         Min.   :0   Length:236         Length:236        
##  Class :character   1st Qu.:0   Class :character   Class :character  
##  Mode  :character   Median :0   Mode  :character   Mode  :character  
##                     Mean   :0                                        
##                     3rd Qu.:0                                        
##                     Max.   :0                                        
##   ENGTYPE_2          NL_NAME_2          VARNAME_2           Shape_Leng    
##  Length:236         Length:236         Length:236         Min.   :0.2102  
##  Class :character   Class :character   Class :character   1st Qu.:0.4752  
##  Mode  :character   Mode  :character   Mode  :character   Median :0.6094  
##                                                           Mean   :0.6897  
##                                                           3rd Qu.:0.7957  
##                                                           Max.   :2.4245  
##    Shape_Area           FID_1           ADM3CODE         ADM1CODE    
##  Min.   :0.001370   Min.   :  0.00   Min.   :  1.00   Min.   :1.000  
##  1st Qu.:0.006363   1st Qu.: 60.75   1st Qu.: 61.75   1st Qu.:2.000  
##  Median :0.010850   Median :196.50   Median :197.50   Median :6.000  
##  Mean   :0.015412   Mean   :161.38   Mean   :162.38   Mean   :4.737  
##  3rd Qu.:0.018074   3rd Qu.:261.25   3rd Qu.:262.25   3rd Qu.:7.000  
##  Max.   :0.087908   Max.   :322.00   Max.   :323.00   Max.   :9.000  
##    ADM1NAME           ADM2NAME            ADM2CODE      ADM3NAME        
##  Length:236         Length:236         Min.   :1103   Length:236        
##  Class :character   Class :character   1st Qu.:2208   Class :character  
##  Mode  :character   Mode  :character   Median :6110   Mode  :character  
##                                        Mean   :4922                     
##                                        3rd Qu.:7175                     
##                                        Max.   :9233                     
##      NUMHH           NUMPER          NUMPER_U         NUMPER_R     
##  Min.   : 2990   Min.   : 12450   Min.   :     0   Min.   :     0  
##  1st Qu.: 8242   1st Qu.: 32880   1st Qu.:     0   1st Qu.: 29711  
##  Median :12660   Median : 51391   Median :     0   Median : 45347  
##  Mean   :15745   Mean   : 65833   Mean   :  9299   Mean   : 52848  
##  3rd Qu.:18097   3rd Qu.: 74626   3rd Qu.:     0   3rd Qu.: 61871  
##  Max.   :86581   Max.   :377396   Max.   :377396   Max.   :209741  
##    NUMPER_EST          NUMPER_M         NUMPER_F         NUMPER17     
##  Min.   :     0.0   Min.   :  6508   Min.   :  5942   Min.   :  4184  
##  1st Qu.:     0.0   1st Qu.: 16170   1st Qu.: 16622   1st Qu.: 11800  
##  Median :    92.5   Median : 25566   Median : 26697   Median : 17878  
##  Mean   :  3685.6   Mean   : 32542   Mean   : 33291   Mean   : 21475  
##  3rd Qu.:  2173.0   3rd Qu.: 36353   3rd Qu.: 37754   3rd Qu.: 24637  
##  Max.   :142483.0   Max.   :203606   Max.   :173790   Max.   :113606  
##     ETSINPOP         ETSLTPOP           ETHINTPO         ETHSLMPO     
##  Min.   :  9026   Min.   :     0.0   Min.   :     0   Min.   :     0  
##  1st Qu.: 28267   1st Qu.:   137.8   1st Qu.:    15   1st Qu.:    76  
##  Median : 44298   Median :   681.0   Median :   284   Median :  1344  
##  Mean   : 54961   Mean   :  2477.9   Mean   :  3484   Mean   :  4464  
##  3rd Qu.: 63581   3rd Qu.:  1679.8   3rd Qu.:  2070   3rd Qu.:  4850  
##  Max.   :206081   Max.   :117510.0   Max.   :141066   Max.   :120109  
##     ETHBURPO          ETHMALPO          ETHOTHPO         RELBUDPO     
##  Min.   :   0.00   Min.   :   0.00   Min.   :   0.0   Min.   :  4463  
##  1st Qu.:   3.00   1st Qu.:   1.75   1st Qu.:   2.0   1st Qu.: 27427  
##  Median :  11.50   Median :   9.00   Median :  11.5   Median : 42501  
##  Mean   : 137.16   Mean   : 192.06   Mean   : 117.4   Mean   : 51381  
##  3rd Qu.:  49.75   3rd Qu.:  47.00   3rd Qu.:  50.5   3rd Qu.: 61387  
##  Max.   :3380.00   Max.   :8680.00   Max.   :5667.0   Max.   :195882  
##     RELHINPO         RELISLPO           RELCATPO          RELOTHPO      
##  Min.   :     0   Min.   :     1.0   Min.   :   11.0   Min.   :   0.00  
##  1st Qu.:   119   1st Qu.:   115.5   1st Qu.:  211.5   1st Qu.:   1.00  
##  Median :   950   Median :  1432.0   Median :  620.5   Median :   7.50  
##  Mean   :  4966   Mean   :  4764.7   Mean   : 4685.6   Mean   :  35.06  
##  3rd Qu.:  3735   3rd Qu.:  5223.0   3rd Qu.: 2064.5   3rd Qu.:  25.00  
##  Max.   :140176   Max.   :134271.0   Max.   :98254.0   Max.   :1217.00  
##     NOHOSU20        NOCLQ200         NOINS200         NONHU200    
##  Min.   : 3651   Min.   :   0.0   Min.   : 14.00   Min.   :  352  
##  1st Qu.: 9311   1st Qu.:  13.0   1st Qu.: 47.00   1st Qu.: 1123  
##  Median :14185   Median :  28.5   Median : 70.50   Median : 1763  
##  Mean   :17071   Mean   : 144.0   Mean   : 85.57   Mean   : 2276  
##  3rd Qu.:19943   3rd Qu.:  67.5   3rd Qu.:109.25   3rd Qu.: 2610  
##  Max.   :66189   Max.   :9029.0   Max.   :628.00   Max.   :27385  
##     TOTBN200        AGHOLB40        AGHLDA40        TOTAGHOL    
##  Min.   : 4173   Min.   :    0   Min.   :    0   Min.   :    0  
##  1st Qu.:10651   1st Qu.: 1087   1st Qu.: 4240   1st Qu.: 7054  
##  Median :15949   Median : 2976   Median : 6004   Median :10300  
##  Mean   :19577   Mean   : 4940   Mean   : 6500   Mean   :11440  
##  3rd Qu.:22854   3rd Qu.: 6127   3rd Qu.: 8441   3rd Qu.:13769  
##  Max.   :95311   Max.   :32141   Max.   :22477   Max.   :42982  
##     NOPNOAL2         NOPOOHG2       NOPOHGOL        NOPOOLO2      
##  Min.   :   0.0   Min.   :   0   Min.   :    0   Min.   :    0.0  
##  1st Qu.: 152.2   1st Qu.: 687   1st Qu.: 1002   1st Qu.:  472.5  
##  Median : 312.0   Median :1248   Median : 1888   Median : 1286.5  
##  Mean   : 468.1   Mean   :1505   Mean   : 2247   Mean   : 2126.9  
##  3rd Qu.: 612.2   3rd Qu.:2146   3rd Qu.: 3169   3rd Qu.: 3174.8  
##  Max.   :5443.0   Max.   :8418   Max.   :10052   Max.   :11436.0  
##     NOPTOT20        NAHLT1AC       AAHLT1AC         NAHB1A2A       ACAHB1A2   
##  Min.   :    0   Min.   :   0   Min.   :   0.0   Min.   :   0   Min.   :   0  
##  1st Qu.: 4173   1st Qu.:1111   1st Qu.: 598.0   1st Qu.:1184   1st Qu.:1524  
##  Median : 6016   Median :1884   Median : 982.5   Median :1810   Median :2393  
##  Mean   : 6347   Mean   :2144   Mean   :1116.7   Mean   :2007   Mean   :2605  
##  3rd Qu.: 8168   3rd Qu.:2926   3rd Qu.:1578.5   3rd Qu.:2738   3rd Qu.:3536  
##  Max.   :20992   Max.   :6787   Max.   :3626.0   Max.   :7033   Max.   :9077  
##     NAHA2AC2        ACAHA2AC        TNAH2001        TACAH200    
##  Min.   :    0   Min.   :    0   Min.   :    0   Min.   :    0  
##  1st Qu.: 1051   1st Qu.: 4230   1st Qu.: 4220   1st Qu.: 7242  
##  Median : 1808   Median : 7221   Median : 6156   Median :11898  
##  Mean   : 2328   Mean   : 8734   Mean   : 6472   Mean   :12456  
##  3rd Qu.: 3222   3rd Qu.:11911   3rd Qu.: 8408   3rd Qu.:17118  
##  Max.   :10840   Max.   :36941   Max.   :21065   Max.   :46643  
##     TEAAC200          RUBAC200         COCACIP2          COCACNIP      
##  Min.   :    0.0   Min.   :   0.0   Min.   :    0.0   Min.   :   0.00  
##  1st Qu.:    0.0   1st Qu.:   0.0   1st Qu.:  191.2   1st Qu.:  40.75  
##  Median :    0.0   Median :   7.0   Median :  510.0   Median : 153.00  
##  Mean   :  929.9   Mean   : 488.5   Mean   : 1628.2   Mean   : 306.21  
##  3rd Qu.:  694.8   3rd Qu.: 245.5   3rd Qu.: 1476.2   3rd Qu.: 437.75  
##  Max.   :11253.0   Max.   :7262.0   Max.   :17939.0   Max.   :2599.00  
##     CASHAC20         TNONCA20          TNOCIM20         DCMPROD2      
##  Min.   :   0.0   Min.   :    0.0   Min.   :   0.0   Min.   :    0.0  
##  1st Qu.:   2.0   1st Qu.:  514.2   1st Qu.: 114.5   1st Qu.:  216.0  
##  Median :  20.5   Median : 1276.5   Median : 237.0   Median :  571.5  
##  Mean   : 155.9   Mean   : 2714.2   Mean   : 446.3   Mean   : 1003.2  
##  3rd Qu.: 177.8   3rd Qu.: 4022.8   3rd Qu.: 587.0   3rd Qu.: 1469.0  
##  Max.   :4369.0   Max.   :20028.0   Max.   :2902.0   Max.   :11254.0  
##     TNOB2001         TNOBIM20         DBMPROD2         TNOGS200      
##  Min.   :   0.0   Min.   :   0.0   Min.   :   0.0   Min.   :   0.00  
##  1st Qu.: 163.8   1st Qu.:   4.0   1st Qu.:   7.0   1st Qu.:  87.75  
##  Median : 385.0   Median :  40.5   Median :  88.0   Median : 229.00  
##  Mean   : 923.3   Mean   : 154.5   Mean   : 253.2   Mean   : 526.47  
##  3rd Qu.:1180.2   3rd Qu.: 146.0   3rd Qu.: 307.8   3rd Qu.: 663.00  
##  Max.   :6389.0   Max.   :1675.0   Max.   :2730.0   Max.   :3654.00  
##     TNOS2001          TNOC2001         SAM1000R         SAM700R2    
##  Min.   :   0.00   Min.   :     0   Min.   :  0.00   Min.   :  180  
##  1st Qu.:   3.75   1st Qu.:  3347   1st Qu.:  4.00   1st Qu.: 1490  
##  Median :  44.50   Median :  7886   Median : 17.00   Median : 2864  
##  Mean   : 229.68   Mean   : 32453   Mean   : 34.59   Mean   : 3320  
##  3rd Qu.: 234.25   3rd Qu.: 26224   3rd Qu.: 44.25   3rd Qu.: 4588  
##  Max.   :3631.00   Max.   :606821   Max.   :379.00   Max.   :11675  
##     SAM500R2         SAM400R2       SAM350R2       SAM250R2     
##  Min.   :  0.00   Min.   :   0   Min.   : 202   Min.   : 111.0  
##  1st Qu.:  2.00   1st Qu.:   3   1st Qu.: 695   1st Qu.: 454.0  
##  Median : 11.00   Median :  13   Median :1080   Median : 772.5  
##  Mean   : 20.29   Mean   :1275   Mean   :1404   Mean   : 889.3  
##  3rd Qu.: 29.00   3rd Qu.:2534   3rd Qu.:1682   3rd Qu.:1157.5  
##  Max.   :220.00   Max.   :9900   Max.   :7077   Max.   :4287.0  
##     SAM140R2          SAMTNOF2        NOPAPAR3        ASSWMAJ3    
##  Min.   :   0.00   Min.   : 1635   Min.   :    0   Min.   :    0  
##  1st Qu.:   0.00   1st Qu.: 4578   1st Qu.: 2528   1st Qu.:    0  
##  Median :   0.00   Median : 6304   Median : 4508   Median :    0  
##  Mean   :  31.03   Mean   : 6961   Mean   : 5422   Mean   : 1948  
##  3rd Qu.:  14.25   3rd Qu.: 8683   3rd Qu.: 7922   3rd Qu.: 1413  
##  Max.   :1535.00   Max.   :19405   Max.   :21754   Max.   :31780  
##     ASSWMIN3          ASSWRF3          ASSWTOT3        PHAMJI00      
##  Min.   :    0.0   Min.   :   0.0   Min.   :    0   Min.   :    0.0  
##  1st Qu.:  284.8   1st Qu.: 140.5   1st Qu.: 1719   1st Qu.:    0.0  
##  Median :  780.0   Median : 746.5   Median : 3256   Median :    0.0  
##  Mean   : 1430.0   Mean   :1302.7   Mean   : 4689   Mean   : 1401.5  
##  3rd Qu.: 1763.2   3rd Qu.:2010.8   3rd Qu.: 5613   3rd Qu.:  587.2  
##  Max.   :12545.0   Max.   :6804.0   Max.   :31780   Max.   :29104.0  
##     PHAMII00          PHARF00           PHATO00           PHAMJI01    
##  Min.   :    0.0   Min.   :   0.00   Min.   :    0.0   Min.   :    0  
##  1st Qu.:  106.5   1st Qu.:  29.75   1st Qu.:  926.5   1st Qu.:    0  
##  Median :  427.5   Median : 422.50   Median : 2005.0   Median :    0  
##  Mean   :  878.4   Mean   : 972.22   Mean   : 3252.2   Mean   : 2797  
##  3rd Qu.: 1096.0   3rd Qu.:1623.50   3rd Qu.: 4090.8   3rd Qu.: 1075  
##  Max.   :12469.0   Max.   :5417.00   Max.   :29104.0   Max.   :59047  
##     PHAMII01        PHARF01        PHATO01         DISROADS     
##  Min.   :    0   Min.   :   0   Min.   :    0   Min.   : 0.050  
##  1st Qu.:  216   1st Qu.:  37   1st Qu.: 1559   1st Qu.: 0.620  
##  Median :  828   Median : 577   Median : 3474   Median : 1.050  
##  Mean   : 1419   Mean   :1423   Mean   : 5640   Mean   : 1.521  
##  3rd Qu.: 1919   3rd Qu.:2293   3rd Qu.: 6873   3rd Qu.: 1.900  
##  Max.   :13522   Max.   :8844   Max.   :59047   Max.   :11.940  
##     DISTOWN          MAHARF          YALARF           TOTRF     
##  Min.   : 0.00   Min.   :143.0   Min.   : 727.0   Min.   :1178  
##  1st Qu.:10.00   1st Qu.:376.0   1st Qu.: 973.8   1st Qu.:1739  
##  Median :18.00   Median :558.0   Median :1049.0   Median :1960  
##  Mean   :19.42   Mean   :564.7   Mean   :1073.4   Mean   :2044  
##  3rd Qu.:27.00   3rd Qu.:745.0   3rd Qu.:1195.0   3rd Qu.:2402  
##  Max.   :48.00   Max.   :968.0   Max.   :1405.0   Max.   :2979  
##     AVE_ELEV         MAX_ELEV           AVE_SLOP           MAX_SLOP     
##  Min.   : 0.155   Min.   :   5.878   Min.   :   2.962   Min.   : 0.976  
##  1st Qu.: 1.171   1st Qu.: 126.250   1st Qu.:  37.417   1st Qu.:21.000  
##  Median : 3.397   Median : 365.500   Median : 101.088   Median :35.000  
##  Mean   : 5.561   Mean   : 616.380   Mean   : 230.341   Mean   :39.022  
##  3rd Qu.: 9.343   3rd Qu.: 950.250   3rd Qu.: 287.656   3rd Qu.:53.500  
##  Max.   :18.785   Max.   :2503.000   Max.   :1670.571   Max.   :85.000  
##     PCTPOOHH         POVGROUP        NUMPOHH          FGT_0       
##  Min.   : 1.001   Min.   :1.000   Min.   :  367   Min.   : 1.143  
##  1st Qu.:23.370   1st Qu.:3.000   1st Qu.: 2401   1st Qu.:28.167  
##  Median :29.555   Median :4.000   Median : 3386   Median :34.892  
##  Mean   :27.737   Mean   :3.962   Mean   : 3698   Mean   :32.614  
##  3rd Qu.:32.803   3rd Qu.:5.000   3rd Qu.: 4707   3rd Qu.:38.345  
##  Max.   :45.650   Max.   :6.000   Max.   :10215   Max.   :52.854  
##     POORPER         PCTSAMHH          SAMBENGP        SAMMISAL    
##  Min.   : 1531   Min.   :  5.197   Min.   :1.000   Min.   :1.000  
##  1st Qu.: 9659   1st Qu.: 37.708   1st Qu.:5.000   1st Qu.:1.000  
##  Median :14187   Median : 48.712   Median :6.000   Median :3.000  
##  Mean   :15406   Mean   : 48.371   Mean   :5.475   Mean   :2.581  
##  3rd Qu.:20407   3rd Qu.: 58.604   3rd Qu.:6.000   3rd Qu.:4.000  
##  Max.   :45218   Max.   :100.525   Max.   :6.000   Max.   :5.000  
##     POPDENS            id2           SMHOLDPAG        AGOPWOLA      
##  Min.   :  8.61   Min.   :  1.00   Min.   : 0.00   Min.   :0.00000  
##  1st Qu.: 68.83   1st Qu.: 61.75   1st Qu.:18.48   1st Qu.:0.03330  
##  Median :131.28   Median :197.50   Median :37.65   Median :0.05320  
##  Mean   :168.03   Mean   :162.38   Mean   :35.50   Mean   :0.07213  
##  3rd Qu.:210.91   3rd Qu.:262.25   3rd Qu.:48.97   3rd Qu.:0.09084  
##  Max.   :778.22   Max.   :323.00   Max.   :88.55   Max.   :0.57604  
##     PACSMALL         NUMAGHH            POOR                 geometry   POORFac
##  Min.   : 0.000   Min.   :0.0000   Min.   :0.0000   MULTIPOLYGON :236   0:179  
##  1st Qu.: 4.524   1st Qu.:0.3197   1st Qu.:0.0000   epsg:4326    :  0   1: 57  
##  Median :11.544   Median :0.5433   Median :0.0000   +proj=long...:  0          
##  Mean   :11.559   Mean   :0.5379   Mean   :0.2415                              
##  3rd Qu.:16.183   3rd Qu.:0.7564   3rd Qu.:0.0000                              
##  Max.   :42.813   Max.   :1.2079   Max.   :1.0000
library(spatialEco)
## Warning: package 'spatialEco' was built under R version 4.4.2
library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
#run a logistic regression
sl_lmodel <- logistic.regression(sl, y="POOR",
x=c("NUMAGHH","DISTOWN","DISROADS", "ASSWMAJ3"))
sl_lmodel$model
## Logistic Regression Model
## 
## rms::lrm(formula = form, data = ldata, x = TRUE, y = TRUE, penalty = bf$penalty)
## 
## 
## Penalty factors
## 
##  simple nonlinear interaction nonlinear.interaction
##       1         1           1                     1
## 
##                       Model Likelihood      Discrimination    Rank Discrim.    
##                             Ratio Test             Indexes          Indexes    
## Obs           236     LR chi2    32.74      R2       0.189    C       0.744    
##  0            179     d.f.       3.851      R2(4,236)0.112    Dxy     0.488    
##  1             57    Pr(> chi2)<0.0001    R2(4,129.7)0.194    gamma   0.488    
## max |deriv| 6e-11     Penalty     0.81      Brier    0.160    tau-a   0.179    
## 
##           Coef    S.E.   Wald Z Pr(>|Z|) Penalty Scale
## Intercept -3.1991 0.4818 -6.64  <0.0001     0.0000    
## NUMAGHH    3.0747 0.8022  3.83  0.0001      0.2773    
## DISTOWN    0.0077 0.0176  0.44  0.6607     11.0418    
## DISROADS   0.0988 0.1149  0.86  0.3899      1.5223    
## ASSWMAJ3  -0.0001 0.0000 -1.27  0.2038   4331.3029
sl_lmodel$diagTable
##          Names        Value
## 1          Obs 2.360000e+02
## 2    Max Deriv 6.355094e-11
## 3   Model L.R. 3.273924e+01
## 4         d.f. 3.851091e+00
## 5            P 1.124815e-06
## 6            C 7.438008e-01
## 7          Dxy 4.876017e-01
## 8        Gamma 4.876017e-01
## 9        Tau-a 1.794086e-01
## 10          R2 1.891640e-01
## 11     R2(236) 1.265524e-01
## 12   R2(4,236) 1.116221e-01
## 13   R2(129.7) 2.182375e-01
## 14 R2(4,129.7) 1.937519e-01
## 15       Brier 1.598022e-01
## 16           g 1.061064e+00
## 17          gr 2.889443e+00
## 18          gp 1.731764e-01
## 19         PEN 1.000000e+00
## 20         AIC 2.503706e+01
sl_lmodel$coefTable
##    Variable          Coef     StdError       Wald         Prob
## 1 Intercept -3.199070e+00 4.817686e-01 -6.6402618 3.131265e-11
## 2   NUMAGHH  3.074701e+00 8.022077e-01  3.8327994 1.266933e-04
## 3   DISTOWN  7.747269e-03 1.764783e-02  0.4389926 6.606669e-01
## 4  DISROADS  9.879432e-02 1.149143e-01  0.8597215 3.899426e-01
## 5  ASSWMAJ3 -5.152669e-05 4.054994e-05 -1.2706970 2.038365e-01
#save residuals and predicted probabilities
VarNamesLogModel <- rownames(sl_lmodel$model$var)[-1]
sl$LogResiduals <- sl_lmodel$Residuals[,1]
sl$LogStdResiduals <- sl_lmodel$Residuals[,2]
sl$LogProbs <- predict(sl_lmodel$model,
 sf::st_drop_geometry(sl[,VarNamesLogModel]),
type='fitted')
#plot the probability of being poor predicted by autologistic model
plot(sl["LogProbs"],main='Poverty probability predicted by Secretariat
Logistic Model', breaks="jenks")

#make a dataframe of the variables in the layer
sl_coord <- st_coordinates(st_centroid(sl))
## Warning: st_centroid assumes attributes are constant over geometries
sl_df <- data.frame(sl)
sl_autolmodel <- logistic.regression(sl_df,
 y="POOR",
x=c("NUMAGHH","DISTOWN","DISROADS", "ASSWMAJ3"),
autologistic=TRUE,
coords=sl_coord)
sl_autolmodel$model
## Logistic Regression Model
## 
## rms::lrm(formula = form, data = ldata, x = TRUE, y = TRUE, penalty = bf$penalty)
## 
## 
## Penalty factors
## 
##  simple nonlinear interaction nonlinear.interaction
##     0.9       0.9         0.9                   0.9
## 
##                       Model Likelihood      Discrimination    Rank Discrim.    
##                             Ratio Test             Indexes          Indexes    
## Obs           236     LR chi2   121.14      R2       0.587    C       0.929    
##  0            179     d.f.       4.698      R2(5,236)0.379    Dxy     0.859    
##  1             57    Pr(> chi2)<0.0001    R2(5,129.7)0.580    gamma   0.859    
## max |deriv| 3e-10     Penalty     3.52      Brier    0.090    tau-a   0.316    
## 
##           Coef    S.E.   Wald Z Pr(>|Z|) Penalty Scale
## Intercept -4.8686 0.7517 -6.48  <0.0001     0.0000    
## NUMAGHH    3.8686 1.1322  3.42  0.0006      0.2631    
## DISTOWN   -0.0138 0.0243 -0.57  0.5699     10.4752    
## DISROADS   0.0997 0.1622  0.61  0.5387      1.4442    
## ASSWMAJ3  -0.0001 0.0001 -0.95  0.3428   4109.0347    
## AutoCov    4.9305 0.6500  7.59  <0.0001     0.3142
sl_autolmodel$diagTable
##          Names        Value
## 1          Obs 2.360000e+02
## 2    Max Deriv 3.224159e-10
## 3   Model L.R. 1.211375e+02
## 4         d.f. 4.698289e+00
## 5            P 0.000000e+00
## 6            C 9.294325e-01
## 7          Dxy 8.588650e-01
## 8        Gamma 8.588650e-01
## 9        Tau-a 3.160115e-01
## 10          R2 5.866516e-01
## 11     R2(236) 3.924751e-01
## 12   R2(5,236) 3.794665e-01
## 13   R2(129.7) 5.961928e-01
## 14 R2(5,129.7) 5.803217e-01
## 15       Brier 8.995618e-02
## 16           g 2.386243e+00
## 17          gr 1.087257e+01
## 18          gp 2.998143e-01
## 19         PEN 9.000000e-01
## 20         AIC 1.117409e+02
sl_autolmodel$coefTable
##    Variable          Coef     StdError       Wald         Prob
## 1 Intercept -4.868592e+00 7.516760e-01 -6.4769823 9.357515e-11
## 2   NUMAGHH  3.868608e+00 1.132242e+00  3.4167693 6.336895e-04
## 3   DISTOWN -1.382178e-02 2.432633e-02 -0.5681819 5.699115e-01
## 4  DISROADS  9.972160e-02 1.622130e-01  0.6147570 5.387152e-01
## 5  ASSWMAJ3 -5.237384e-05 5.520869e-05 -0.9486520 3.427976e-01
## 6   AutoCov  4.930545e+00 6.499581e-01  7.5859426 3.300784e-14
#add the residual etc to the dataframe
sl$AutoCov <- sl_autolmodel$AutoCov
sl$AutoLogResiduals <- sl_autolmodel$Residuals[,1]
sl$AutoLogStdResiduals <- sl_autolmodel$Residuals[,2]
sl$AutoLogProbs <- predict(sl_autolmodel$model,
 type='fitted')
#plot the probability of being poor predicted by autologistic model
plot(sl["AutoLogProbs"],main='Poverty probability predicted by Secretariat
Autologistic Model', breaks="jenks")

OLS and GWR in R

library(spdep) 
library(spgwr) 
## Warning: package 'spgwr' was built under R version 4.4.3
## Loading required package: sp
## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)
library(RColorBrewer)
library(tmap)
## Warning: package 'tmap' was built under R version 4.4.3
# OLS model
sl_pov_olsmodel <- lm(PCTPOOHH ~ NUMAGHH + MAHARF + DISROADS + DISTOWN + ASSWMAJ3, data=sl) 
summary(sl_pov_olsmodel)
## 
## Call:
## lm(formula = PCTPOOHH ~ NUMAGHH + MAHARF + DISROADS + DISTOWN + 
##     ASSWMAJ3, data = sl)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.5146  -4.3922  -0.1625   3.1168  17.1759 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 24.6716175  1.9800545  12.460  < 2e-16 ***
## NUMAGHH     17.5793199  1.9896718   8.835 2.62e-16 ***
## MAHARF      -0.0111072  0.0024187  -4.592 7.22e-06 ***
## DISROADS    -0.2471006  0.3511222  -0.704   0.4823    
## DISTOWN      0.0389303  0.0455097   0.855   0.3932    
## ASSWMAJ3    -0.0002558  0.0001085  -2.358   0.0192 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.142 on 230 degrees of freedom
## Multiple R-squared:  0.487,  Adjusted R-squared:  0.4758 
## F-statistic: 43.66 on 5 and 230 DF,  p-value: < 2.2e-16
AIC(sl_pov_olsmodel)  
## [1] 1534.377
# GWR models 
adaptive_bw <- gwr.sel(PCTPOOHH ~ NUMAGHH + MAHARF + DISROADS + 
                         DISTOWN + ASSWMAJ3, 
                       data=sl, method = "aic", adapt = TRUE, coords=sl_coord) 
## Bandwidth: 0.381966 AIC: 1480.366 
## Bandwidth: 0.618034 AIC: 1502.147 
## Bandwidth: 0.236068 AIC: 1450.837 
## Bandwidth: 0.145898 AIC: 1411.916 
## Bandwidth: 0.09016994 AIC: 1370.652 
## Bandwidth: 0.05572809 AIC: 1332.262 
## Bandwidth: 0.03444185 AIC: 1321.457 
## Bandwidth: 0.02178569 AIC: 1347.974 
## Bandwidth: 0.04177523 AIC: 1318.577 
## Bandwidth: 0.04115187 AIC: 1319 
## Bandwidth: 0.04710475 AIC: 1321.653 
## Bandwidth: 0.04381093 AIC: 1318.699 
## Bandwidth: 0.04268552 AIC: 1318.296 
## Bandwidth: 0.04272621 AIC: 1318.304 
## Bandwidth: 0.04252023 AIC: 1318.269 
## Bandwidth: 0.04223567 AIC: 1318.319 
## Bandwidth: 0.04241153 AIC: 1318.255 
## Bandwidth: 0.04234436 AIC: 1318.265 
## Bandwidth: 0.04245222 AIC: 1318.26 
## Bandwidth: 0.04241153 AIC: 1318.255
adaptive_bw 
## [1] 0.04241153
# Optimal adaptive bandwidth for the above case is: approximately 0.04
# each regression will expand or contract to include roughly 4% of sample

fixed_bw <- ggwr.sel(PCTPOOHH ~ NUMAGHH + MAHARF + DISROADS + 
                       DISTOWN + ASSWMAJ3, 
                     data=sl, adapt = FALSE, coords=sl_coord)
## Bandwidth: 1.297489 CV score: 8242.629 
## Bandwidth: 2.097285 CV score: 8755.026 
## Bandwidth: 0.8031877 CV score: 7255.535 
## Bandwidth: 0.4976927 CV score: 5854.028 
## Bandwidth: 0.3088864 CV score: 4421.947 
## Bandwidth: 0.1921977 CV score: 3786.395 
## Bandwidth: 0.1200801 CV score: 4290.592 
## Bandwidth: 0.2092028 CV score: 3803.537 
## Bandwidth: 0.1950849 CV score: 3786.397 
## Bandwidth: 0.1936372 CV score: 3786.24 
## Bandwidth: 0.1936779 CV score: 3786.24 
## Bandwidth: 0.1935965 CV score: 3786.24 
## Bandwidth: 0.1936372 CV score: 3786.24
fixed_bw
## [1] 0.1936372
# 21415.7 
# This means that the bandwidth for each regression using a fixed bandwidth 
# will include all centroid within approx 21 km  


# USING THE ADAPTIVE BISQUARE KERNEL 
# Note that the input for the adapt command is the value generated 
# in the calculation above for adaptive.bw  
sl_adgwr_model <- gwr(PCTPOOHH ~ NUMAGHH + MAHARF + DISROADS + DISTOWN + 
                     ASSWMAJ3, data=sl, adapt = adaptive_bw, 
                     hatmatrix = TRUE,
                     coords=sl_coord) 
names(sl_adgwr_model) 
##  [1] "SDF"       "lhat"      "lm"        "results"   "bandwidth" "adapt"    
##  [7] "hatmatrix" "gweight"   "gTSS"      "this.call" "fp.given"  "timings"
summary(sl_adgwr_model) 
##           Length Class                  Mode     
## SDF         236  SpatialPointsDataFrame S4       
## lhat      55696  -none-                 numeric  
## lm           11  -none-                 list     
## results      14  -none-                 list     
## bandwidth   236  -none-                 numeric  
## adapt         1  -none-                 numeric  
## hatmatrix     1  -none-                 logical  
## gweight       1  -none-                 character
## gTSS          1  -none-                 numeric  
## this.call     6  -none-                 call     
## fp.given      1  -none-                 logical  
## timings      12  -none-                 numeric
names(sl_adgwr_model$results) 
##  [1] "v1"       "v2"       "delta1"   "delta2"   "sigma2"   "sigma2.b"
##  [7] "AICb"     "AICh"     "AICc"     "edf"      "rss"      "nu1"     
## [13] "odelta2"  "n"
sl_adgwr_model$results$AICc
## [1] 1330.992
#Note the drop in AICc compared to regression AIC

# save local coefficients and standard errors to a new dataframe
gwrmodel_results <- as.data.frame(sl_adgwr_model$SDF)

# join it with the shapefile using ID_2 to match rows for both
# join it with the spatial data 
sl_gwr <- cbind(sl, as.matrix(gwrmodel_results))
sl_gwr_sf <- st_as_sf(sl_gwr)


# write your spatial data from R to a shapefile 
write_sf(sl_gwr, "SriLankaGWR.shp", OVERWRITE=TRUE)
## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver
# Repeat this section for each variable in your GWR
# first create a new variable NUMAGHH_T with significant coeff
#note that R has added a .1 after the variable to show that 
#this is the coefficient and not the original variable
sl_gwr$NUMAGHH_T <- sl_gwr$NUMAGHH.1/sl_gwr$NUMAGHH_se
sl_gwr$DISROADS_T <- sl_gwr$DISROADS.1/sl_gwr$DISROADS_se

#If there are values of Max and/or Min over 1.96 we map it
#if its not significant use the OLS model coefficient 
summary(sl_gwr$NUMAGHH_T)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.6777  2.3531  3.3901  4.0401  5.3722 13.5991
summary(sl_gwr$DISROADS_T)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -3.5373 -0.0109  0.4778  0.6146  1.4391  5.9105
#setting the coefficient so that it is only symbolized if T is significant
sl_gwr$CoeffNUMAGHH <- ifelse(abs(sl_gwr$NUMAGHH_T)>=1.96,sl_gwr$NUMAGHH_T, NA)
sl_gwr$CoeffDISROADS <- ifelse(abs(sl_gwr$DISROADS_T)>=1.96,sl_gwr$DISROADS_T, NA)


#Tmap is better for mapping multiple layers
#tmap_mode is set to view so that you can zoom in to see location names 
#change the basemap from the default to Openstreetmap to see place names

#see all pallettes by running the next line without the comment #
#cols4all::c4a_palettes() 


map_slr2 <- tm_shape(sl_gwr) +
                tm_polygons(fill="localR2", 
                            fill.scale = tm_scale_intervals(values="brewer.blues", style = "jenks"),
                            fill.legend = tm_legend(title = "R Squared", size = 0.8))+ 
                tm_borders() + tmap_mode("view")
## ℹ tmap mode set to "view".
map_slr2